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Abstract 



1 Introduction 



An extension of the SARGE-algorithm of Q is introduced, which includes the in- 

^ ' coming momenta in the kinematical pole structure of the density with which the 

■ momenta are generated. The algorithm is compared with RAMBO in the integration 
' of QCD-amplitudes in the SPHEL-approximation, and the computing times are ex- 
. trapolated to those for the calculation with exact matrix elements. 

o 
o 

:a 

Qh! In future experiments with hadron colliders, such as the LHC, many multi-jet events will 
r-| ! occur. These can be divided into interesting events (IE), and the background. The main 
difference between the two classes is that the Standard Model shall not have proven yet its 

■ capability of dealing with the description of the IE at the moment when they are analyzed. 
^ I The background shall not manifest itself as such a heavy test for the standard model. 

However, we still need to know the cross sections of the background in order to compare 
the ratio of these and those of the IE with the predictions of the Standard Model. 

For large part of the background, a piece of the transition amplitude consists of a multi- 
parton QCD-amplitude, and it is well known |^ that it contributes to the cross section 
with a singular behavior in phase space (PS), given by the so-called antenna pole structure 
(APS). In particular, for processes involving only n gluons the most important contribution 
comes from the sum of all permutations in the momenta of 

(1) 



{PrP2){P2-P3){P3-P4) ■ ■ ■ iPn-rPn)iPn-Pl 



and the singular nature stems from the fact that the scalar products Pi-pj can become 
very small. For the calculation of integrals over PS, the Monte Carlo (MC) method is the 
only option, so that an algorithm to generate random momenta is needed. For processes at 
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single antenna integrated to 1% error 


number of 


cut-off 


number of 


momenta 


CM-energy 


PS points 


3 


0.183 


10,069 


4 


0.129 


26,401 


5 


0.100 


58, 799 


6 


0.0816 


130,591 


7 


0.0690 


240, 436 


8 


0.0598 


610,570 



evaluation amplitude in 1 PS point 


number of 
final gluons 


cpu-time (seconds) 


SPHEL 


exact 


3 


2.83x10^5 


1.60x10-1 


4 


9.76x10-5 


5.54x10-1 


5 


4.88x10"^ 


1.945 


6 


3.26x10-3 


6.06 


7 


2.57x10-2 


19.91 


8 




64.45 



Table 1: Typical number of PS points and computing times. 



high energy, the momenta may be massless, and RAMBO generates any number of them 
distributed uniformly in PS. This uniform distribution, however, has the disadvantage that, 
for the integration of an integrand containing the APS, a large number of events is needed 
to reach a result to acceptable precision. As an illustration, we present in the left table 
of Tab. |I| the number of PS points needed to integrate the single antenna of Eq. (|T]), so 
not even the sum of its permutations, to an expected error of 1%. The antenna cannot be 
integrated over the whole of PS because of the singularities, so these have to be cut out. 
This is done through the restriction {pi + pj)"^ > sq for all i, j = 1, . . . ,n, and in the table 
the ratio between ^/sq and the total energy ^/s is given. These numbers are based on the 
reasonable choice so/s = 0.2/[n(?T, — 1)]. 

Performing MC integration with very many events is not a problem if the evaluation 
of the integrand in each PS point is cheap in computing time. This is, for example, the 
case for algorithms to calculate the squared multi-parton amplitudes based on the so called 
SPHEL-approximation, for which only the kinematical structure of (|I|) is implemented . 
Nowadays, algorithms to calculate the exact matrix elements exist, which are far more 
time-consuming [0, §]. As an illustration of what is meant by 'more time-consuming', we 
present the right table of Tab.|l| with the typical cpu-time needed for the evaluation in one 
PS point of the integrand for processes of two gluons going to more gluons, both for the 
SPHEL-approximation and the exact matrix elements [|lT|. It is expected, and observed, 
that the exact matrix elements reveal the same kind of singularity structures as the APS, 
so that, according to the tables, the PS integration for a process with 8 final gluons would 
take in the order of 400 days . . . 

The solution to this problem is importance sampling. Instead of RAMBO, a PS generator 
should be used which generates momenta with a density including the APS. In we 
introduced an algorithm that does part of the job, and is called SARGE (from Staggered 
Antenna Radiation Generator). It generates n random momenta with a density proportional 
to the the sum of all permutations of (|l]), and because they are all random, they should 
be interpreted as outgoing momenta. In the pole structure of a "real" QCD-amplitude, 
however, also the incoming momenta occur. In this paper, we introduce an extension of 



2 



the SARGE-algorithm, which includes these pole structures. We compare it with RAMBO in 
the calculation of integrals of QCD-amplitudes in the SPHEL-approximation, and extrap- 
olate the computing times to those for the calculation with exact matrix elements. The 
conclusion will be that SARGE takes account for a substantial reduction in computing time. 

For the sake of completeness, we describe the full algorithm in this paper, including 
the piece that was introduced in 0]. What we actually presented there was only the 
algorithm and no proof of its correctness whatsoever. In this paper, however, we shall 
meet our engagements with respect to this. We do this with the help of the unitary 
algorithm formalism, which we introduce in the following section by its application to the 
RAMBO-algorithm. 

2 Notation and the unitary algorithm formahsm 

The relativistic momentum p = p^jp"^, p^) = p) of an elementary particle is a vector 
in R^. The momentum with the opposite 3- momentum is denoted by 

p = -p) . (2) 

We shall need the first and the fourth canonical basis vectors, which we denote 

Co =(1,0,0,0) and eg = (0, 0, 0, 1) . (3) 

A typical parameterization of a 3-momentum with unit length is given by n{z, ip), where 

ni{z,(p) = Vl — z'^ simp , n2{z, (p) = y/l — z'^ cos , n3{z,ip) = z . (4) 

The Lorentz invariant scalar product shall be denoted with a dot or with parentheses: 

(pq) = p-q = p^q^ — p-q , p-q = p^q^ + p^q^ + p'^q^ . (5) 

The product of a vector with itself is denoted as a square 

p' = {pp) = {p'f-\p? , \P\ = {P-Pf" ■ (6) 

The same notation for the quadratic form and the 2-component will not lead to confusion, 
because the 2-component will not appear explicitly anymore after this section. For physical 
particles, p^ has to be positive, and in that case, the square root gives the invariant mass 
of the particle: 

"^P = if > . (7) 

A boost that transforms a momentum p, with p^ > 0, to mpCo is denoted Tip, so 

TipP = rripeo and mpHpeo = p . (8) 
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A rotation that transforms p to p^eo + is denoted TZp, so 

Kpp = p°eo + |p|e3 and Upp = p^eo - \p\e3 . (9) 

Since rotations only change the 3-momentum, we shall use the same symbol if a rotation 
is restricted to three-dimensional space. 

The physical PS of n particles is the (3n — 4)-dimensional subspace of R"'", given by 
the restrictions that the energies of the particles are positive, the invariant masses squared 
pj are fixed to given positive values Sj, and that the sum 

n 



Pin) ^^Pi (10) 



i=l 



of the momenta is fixed to a given momentum P. The restrictions for the separate momenta 
shall be expressed with a 'PS characteristic distribution' 

^.,(p) = <^(p'-sO^(/) , and m-Mp) ■ (11) 

The generic PS integral, of a function F of a set {p}n = {Pi, • • • )Pn} of momenta, that 
has to be calculated is then given by 

An integral shall always start with a single J -symbol, and for every integration variable, 
say a dz means 'integrate z over the appropriate integration region'. If it is not evident 
what this region is, it shall be made explicit with the help the of logical ^-functions, which 
have statements H as argument, and are defined through 

, ^ I 1 if n is true , , 

^ n = <^ 13 

if n is false. 

2.1 The RAMBO algorithm in the UAF 

RAMBO was developed with the aim to generate the fiat PS distribution of n massless 
momenta as uniformly as possible, and such that the sum of the momenta is equal to 
y^Co with s a given squared energy. This means that the system of momenta is in its 
center-of-mass frame (CMF), and that the density is proportional to the 'PS characteristic 
distribution' 

n 

e.(Mn) = 5\p^n)-^^seo)\[^{p^) ■ (M) 

1=1 

The algorithm consists of the following steps: 
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Algorithm 1 (RAMBO) 

1. generate n massless vectors qj with positive energy without constraints but under 
some normahzed density /(^j); 

2. compute the sum q(^n) of the momenta g^; 

3. determine the Lorentz boost and scahng transform that bring g(„) to a/s^q; 

4. perform these transformations on the gj, and call the result pj. 

Trivially, the algorithm generates momenta that satisfy the various (5-constraints, but it 
is not clear a priori that the momenta have the correct distribution. To prove that they 
actually do, we apply the unitary algorithm formalism (UAF). We write the generation of 
a variable as the integral of the density with which that variable is generated, and interpret 
every assignment as a generation with a density that is given by a Dirac delta-distribution. 
Only the assignment of the final output should not be written as an integral, but only with 
the delta-distributions. The UAF tells us that Algorithm |l| generates a density 



/I L 
A — 1 ^^Q(n') 



j = l 

n 

\[5\p,-xHkqi) . (15) 



X 

1=1 



The unitarity of the algorithm is expressed by the fact that integration of the above equa- 
tion over the set of variables {p}n leads to the identity 1 = 1. To calculate the distribution 
yielded by this algorithm, the integral has to be evaluated. First of all, some simple algebra 
using p(„) = xHbq(n), <l(n) = x^^l-i^^p(n) and the Lorentz and scaling properties of the Dirac 
^-distributions leads to 

5'(h-^)5(x-^) = —5^(p(„)-v^eo) 5(6^-1) . (16) 



Furthermore, since we may write 

1 

X 

under the integral, the l.h.s. of Eq. (|T5]) becomes 



d%5{q])5\p,-xn,q,) = -5{p]) (17) 



/ Qs{{v}n) d%5{h'' - i)dx-^^\{f{-n,'p,)e{e,-n^'p, > o) . 

J X .^^ X 

In the standard RAMBO algorithm, the following choice is made for /: 



f{q) = ^ exp(-cg°) , (19) 
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where c is a positive number with the dimension of an inverse mass. Therefore, if we use 
that = v^eo and that = e^-q for any g, then 

1=1 ^ ^ i=l 



2 \ " 



^ , exp ( -^^6" ) ^(6° > 0) . (20) 

As a result of this, the variables pj, i = 1, . . . ,n only appear in O^, as required. The 
remaining integral is calculated in Appendix A, with the result that RAMBO generates the 
density 

*.({pW^e,(w„)(i)""£(!^H^. pi) 

Incidentally, we have computed here the volume of the PS for n massless particles: 



L 



R4n V2/ r(n)r(n — 1) 

Note, moreover, that c does not appear in the final answer; this is only natural since any 
change in c will automatically be compensated by a change in the computed value for 
X. Finally, it is important to realize that the 'original' PS has dimension 3n, while the 
resulting one has dimension 3n — 4: there are configurations of the momenta qj that are 
different, but after boosting and scaling end up as the same configuration of the pj. It is 
this reduction of the dimensionality that necessitates the integrals over b and x. 



3 The basic antenna 

As mentioned before, we want to generate momenta that represent radiated partons with 
a density that has the antenna structure [(piP2)(p2P3)(p3P4) • • • , ipn-iPn)iPnPi)]''^ ■ Nat- 
urally, the momenta can be viewed as coming from a splitting process: one starts with 
two momenta, a third is radiated off creating a new pair of momenta of which a fourth is 
radiated off and so on. In fact, models similar to this are used in full-fledged Montc-Carlo 
generators like HERWIG. Let us therefore first try to generate a single massless momentum 
k, radiated from a pair of given massless momenta pi and p2. In order for the distribution 
to have the correct infrared and coUinear behavior, it should qualitatively be proportional 
to [{pik){kp2)]~^- Furthermore, we want the density to be invariant under Lorentz trans- 
formations and scaling of the momenta, keeping in mind that the momenta are three out 
of possibly more in a CMF and that we have to perform these transformations in the end, 
like in RAMBO. This motivates us to define the basic antenna structure as 
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Here, g is a function that serves to regularize the infrared and colhnear singularities, as 
well as to ensure normalization over the whole space for k: therefore, g{C,) has to vanish 
sufficiently fast for both ^ ^ and ^ ^ oo. To find out how k could be generated, we 
evaluate / dA in the CMF of pi and p2- Writing 

E = V (piP2)/2 , p = Hp.+p^pi , q = Hp^+p^k , (24) 

we have 

{piP2)^2E^ , {pik) = Eq''{l-z) , {kp2) = Eq^l + z) , (25) 

where z — p-q/{\p\\q\)- The azimuthal angle of q is denoted (f, so that q — \q\TZ~^n{z, </?). 
We can write 



where. 



d%^{k) = U^dq^difdz = UpiP2)d(pd^id^2 , (26) 



6 = /^ and & = i^, (27) 



SO that z = {^2 — Ci)/(C2 + Ci) ^-iid = -£'(^2 + Ci)- The integral over dA takes on the 
particularly simple form 



/ 



dAip,,p2;k) = (^l^^d^^giO^ . (28) 



The antenna dA{pi,p2; k) will therefore correspond to a unitary algorithm when we let the 
density g be normalized by 

l^^^\s{0 = ^- (29) 

Note that the normalization of dA fixes the overall factor uniquely: in particular the 
appearance of the numerator (P1P2) is forced upon us by the unitarity requirement. 

For g we want to take, at this point, the simplest possible function we can think of, 
that has a sufficiently regularizing behavior. We introduce a positive non-zero number 
and take 

^(0 - 7^7^ die' < e < em) . (30) 

2 log 4m 

The number gives a cut-off for the quotients and ^2 of the scalar products of the 
momenta, and not for the scalar products themselves. It is, however, possible to relate 
to the total energy ^/s in the CMF and a cut-off Sq on the invariant masses, i.e., the 
requirement that 

(Pi ^ '^0 for all momenta Pi ^ pj. (31) 
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This can be done by choosing 

. (32) 

With this choice, the invariant masses (pi + fc)^ and ik \ 'P'lf are regularized, but can 
still be smaller than sq so that the whole of PS, cut by (|3T|) , is covered. The Sq can be 
derived from physical cuts on the transverse momenta and 6*0 on the angles between the 
outgoing momenta: 




With this choice, PS with the physical cuts is covered by PS with the cut of (^). To 
generate the physical PS, the method of hit-and-miss Monte Carlo can be used, i.e, if 
momenta of an event do not satisfy the cuts, the whole event is rejected. We end this 
section with the piece of the PS algorithm that corresponds to the basic dA{j)x,'p2\ k): 

Algorithm 2 (BASIC ANTENNA) 

1. given {pi,p2}, put p <— Hp^+p^pi and put E ^ \/TjhP2)J^ ; 

2. generate two numbers ^i, ^2 independently, each from the density g{^)/^, and (p 
uniformly in [0, 27r) ; 

3. put z <- (6 - ei)/(6 + 6), g° <- ^(6 + 6) and q ^ q'TZ^M^, ip) ; 

4. put A; ^ T^pr+P2l ; 



4 A complete QCD antenna 

The straightforward way to generate n momenta with the antenna structured density is by 
repeated use of the basic antenna. Let us denote 

c/^fc = rfA(g,-,gfc;g,) , (34) 

then 
where 

gn{{q}n) = g[-, \]9[i \]9\i \]9\i \]---9[-( \] ■ (35) 

\{qiqn)J \{qiqn)J \[q2qn)J \{q2qn)J \{qn-2qn)J 

So if we have two momenta qi and g„, then we can easily generate n — 2 momenta qj with 
the antenna structure. Remember that this differential PS volume is completely invariant 
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under Lorentz transformations and scaling transformations, so that it seems self-evident 
to force the set of generated momenta in the CMF with a given energy, using the same 
kind of transformation as in the case of RAMBO. If the first two momenta are generated with 
density /(gi, qn), then the UAF tells us that generated density A'^'^°{{p}n) satisfies 

^r°({p}n) = / d\^{q,)d'qn^{qn)f{quqn)dAlJAlJAl^...dA^^Zl^ 

n 

X d%5^{h- q(^r,)/mq^^^)dx5{x - ^/s /rriq^^^) JJ - xT^^gi) . (36) 

i=l 

If we apply the same manipulations as in the proof of the correctness of RAMBO, we obtain 
the equation 



vr" '^{piP2){P2P-i){PzPA) ■ ■ ■ (Pn-lPn) 

X j d% 6{b^ -l)dx'^ f{x-'n^'pi , x-'n^'pn) . (37) 

Now we choose / such that qi and qn are generated back-to-back in their CMF with total 
energy ^/s , i.e., 

2 

f{Qi,Qn) = -S'^iqi + qn- y/seo) . (38) 

TT 

If we evaluate the second line of Eq. (^^ with this /, we arrive at 



TT 



J dx^ d%5{h^ - 1) 6\x-^n^\pi + pn) - v^eo) 



4 fi. ll (39) 

TT Jo X^ y j 27r(pip„)2 

SO that the generated density is given by 

^r(Wn) = 0.(Wn)^7 — ^ — -^^^MA — — ^ . (40) 

27r" 2 {piP2){p2P3){P3P4) " " " [Pn-lPn) {PnPl) 

Note that, somewhat surprisingly, also the factor {pnPi)^^ comes out, thereby making the 
antenna even more symmetric. In fact, if the density f{qi,q2) = c^exp(— cgj* — cg2)/4vr2 
is taken instead of the one we just used, the calculation can again be done exactly, with 
exactly the same result. The algorithm to generate n momenta with the above antenna 
structure is given by 

Algorithm 3 (QCD ANTENNA) 

1. generate massless momenta qi and g„; 

2. generate n — 2 momenta qj by the basic antennas dAl^^dA^^^dA^ ,^ ■ ■ ■ dA^Z^n'^ 
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3. compute q(^n) — YTj=\ Qj^ boost and scaling transforms that bring q(^n) to 
a/sco; 

4. for j' = 1, . . . , n, boost and scale the qj accordingly, into the pj. 

Usually, the event generator is used to generate cut PS. If a generated event does not 
satisfy the physical cuts, it is rejected. In the calculation of the weight coming with an 
event, the only contribution coming from the functions g is, therefore, their normalization. 
In total, this gives a factor 1/(2 log ^m)^""^ in the density. 



5 Incoming momenta and symmetrization 

The density given by the algorithm above, is not quite what we want. First of all, we 
want to include the incoming momenta po and po in the APS, so that the density becomes 
proportional to [(poPi) (pm) ■ ■ ■ iPn-iPn)iPnPo)V^ instead of [{piP2) ■ ■ ■ {Pn~iPn) iPnPi)]~^ ■ 
Then we want the sum of all permutations of the momenta, including the incoming ones. 



5.1 Generating incoming momenta 

The incoming momenta can be generated after the antenna has been generated. To show 
how, let us introduce the following "regularized" scalar product: 

{pq)s = {pq) + 6p\'' , (41) 

where 5 is a small positive number. This regularization is not completely Lorentz invariant, 
but that does not matter here. Important is that it is still invariant under rotations, as 
we shall see. Using this regularization, we are able to generate a momentum k with a 
probability density 

1 m s{k' - 1) ^^2) 



27rIs{pi,P2) {pik)s{kp2)5 

To show how, we calculate the normalization Is{pi,p2)- Using the Feynman- representation 
of l/[{pik)s{kp2)s], it is easy to see that 

If dx 

where p^ = xpi + {x — l)p2- The integral over z and </? can now be performed, with the 
result that 



(l + <)) -IPxI 2(piP2) Jo {x+-x){x-X-) 
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where x± are the solutions for x of the equation 1 + 5 = \pj. Further evaluation finally 
leads to 



h{PuP2) = log 

X-L. X_ 



X+ 




2pM(2*^ (45, 



Notice that there is a smooth limit to the case in which pi and p2 are back-to-back: 

Isip,p) = lim/,(p,g) = (^0)2(2^5 + ^2) ■ (46) 

The algorithm to generate k can be derived by reading the evaluations of the integrals 
backwards. 

Because k and k are back-to-back, they can serve as the incoming momenta. To fix 
them to Co + 63 and eg — 63, the whole system of momenta can be rotated. If we generate 
momenta with the density use the first two momenta to generate the incoming 

momenta and rotate, we get a density 

= ({p}n) i5{Pi,P2r' f d'km s{k' - 1) ^ ^ ' - , (4?) 

J [piTCkk)s{p2Tlkk)s 
where we used the fact that the whole expression is invariant under rotations, and that 
these are orthogonal transformations. The last line of the previous expression can be 
evaluated further with the result that 

Dsi{p}n) - ^r°(Mn) /'^^''f'^ ' With po^eo + ea, po = eo-e3. (48) 

{PlP0)s{P0P2)5 

The algorithm to generate the incoming momenta is given by 

Algorithm 4 (INCOMING MOMENTA) 

1. given a pair {pi,P2}, calculate x^ and 

2. generate x in [0, 1] with density ~ [(x+ — x)(x — x_)]~^, and put p^ <— xpi + {x — l)p2 ; 

3. generate (p uniformly in [0, 27r), z in [—1, 1] with density ~ (1 + 5 — |p^|-2)~^ ; 

4. put k '}Z~^n{z, ifi) and 1 ; 

5. rotate all momenta with TZk ; 

6. put po ^ ^^/s (eo + 63) and po ^ ^^/s (cq - 63) . 

Notice that Is{pi,P2){piPo)s{PoP2)s is invariant under the scaling pi,p2 cpi,cp2 with a 
constant c, so that scaling of po and po has no infiuence on the density. 

The pair {qi, ^2) with which k is generated is free to choose because we want to sym- 
metrize in the end anyway. We should only choose it such, that we get rid of the factor 
(51^2) in the denominator of Af'^°{{q}n). 
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5.2 Choosing the type of antenna with incoming momenta 

A density which is the sum over permutations can be obtained by generating random per- 
mutations, and returning the generated momenta with permutated labels. This, however, 
only makes sense for the outgoing momenta. The incoming momenta are fixed, and should 
be returned separately from the outgoing momenta by the event generator. Therefore, a 
part of the permutations has to be generated explicitly. There are two kinds of terms in 
the sum: those in which (poPo) appears, and those in which it does not. 

Case 1: antenna with (poPo). To generate the first kind, we can choose a label i at 
random with weight where Si({p}„) is the sum of all scalar products in 

the antenna [|: 

n 

Sl({p}n) = J2^P^P^+l) ■ (49) 
i=l 

This is a proper weight, since all scalar products are positive. The total density gets this 
extra factor then, so that (piPi+i) cancels. The denominator of the weight factor does not 
give a problem, because its singular structure is much softer than the one of the antenna. 
The pair {pi,pi+i} can then be used to generate the incoming momenta, as shown above. 
So in this density y4^'^°({p}„)_Bi({p}„)/Si({p}„) is generated, where 

R /J X N {PiPi+l)l5{Pi,Pi + l)'^ 

Bi{{p}n) = > — J ^ • 50 

^ {PiPo)s{PoPi+i)s 

Case 2: antenna without {poPo)- To generate the second kind, we can choose two 
non-equal labels i and j with weight {piPi+i){pjPj+i) /'S2{{p}n) , where 

n 

^2i{p}n) = Y{PiPi+l){PjPj+l) ■ (51) 

Next, a pair {k, I) of labels has to be chosen from the set of pairs 

= {(hj) , {^,J + 1) , (^ + 1, j) , (^ + 1, J + 1)} • (52) 
If this is done with weight Is{pk,Pi)/^i,j{{p}n), where 

^i,j{{p}n) = ^s{Pk,Pi) , (53) 

(fc,/)G{(i,i)} + 

then the density A'^^^ {{p}n)B2{{p}n) /^2{{p}n) is generated, where 

Is{pk,Pi) Is{pk,Pl)~^ 



^i,j{{p}n) iPkPo)s{PoPl)5 



{PiPi+i){PjPj+i) J2(k,i)e{(i,j)}+iPkPo)s{PoPi)5 ^^^^ 

^ {PiPo)5{Pi+lPo)5{poPj)5{poPj+l)5 T.ik,l)em)}+^s{Pk,Pl) 



^Read i + 1 mod n when i + 1 occurs in this section 
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Before all this, we first have to choose between the two cases, and the natural way to do 
this is with relative weights |sEi({j9}„) and S2({p}„), so that the complete density is equal 
to 



perm. ^ ^ \ li j i'- 



^2{{p}n) 



(55) 



where the first sum is over all permutations of (1,... ,n). One can, of course, try to 
optimize the weights for the two cases using the adaptive multichannel method (cf. 0). 
The result of using the sum of the two densities is that the factors [piPi+i) in the numerator 
of Bi{{p}n) and {piPi+i){pjPj+i) in the numerator of B2{{p}n) cancel with the same factors 
in the denominator of A'^'^^{{p}n), so that we get exactly the pole structure we want. The 
'unwanted' singularities in Bi{{p}n), B2{{p}n) and T,2{{p}n) are much softer than 

the ones remaining in A^^^{{p}n), and cause to trouble. The algorithm to generate the 
incoming momenta and the permutation is given by 

Algorithm 5 (CHOOSE INCOMING POLE STRUCTURE) 

1. choose case 1 or 2 with relative weights |sSi({]9}„) and T,2{{p}n) 

2. in case 1, choose ii with relative weight (pj^pj^+i) and put 12 ^ ii + I ; 

3. in case 2, choose (z, j) with {i 7^ j) and relative weight {piPi+i){pjPj+i), and then 
choose {ii,i2) from {(i,j)}+ with relative weight h{phi'Pi2) i 

4. use {pii,Pi2} to generate the incoming momenta with Algorithm Q 

5. generate a random permutation a E Sn and put pi <— Po-(j) for alH = 1, . . . , ra. 

An algorithm to generate the random permutations can be found in [|^. An efficient 
algorithm to calculate a sum over permutations can be found in 0. 



6 Improvements 



When doing calculations with this algorithm on a PS, cut such that [pi + Pjf' > sq for all 
i j and some reasonable Sq > 0, we notice that a very high percentage of the generated 
events does not pass the cuts. An important reason why this happens is that the cuts, 
generated by the choices of g (Eq. (pOD) and C,^ (Eq. (p^)), are implemented only on quotients 
of scalar products that appear explicitly in the generation of the QCD-antenna: 

,-i \Pi—lPi) J \PiPn) no i /crr\ 

= 7 T and ^2 = 7 z = 2,3...,n-l. (56) 

[Pi-lPn) [Pi-lPn) 

The total number of these ^-variables is 

n^ = 2n~A , (57) 
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and the cuts are implemented such that < 2 < for i = 2, 3 . . . , n — 1. We show 
now how these cuts can be implemented on all quotients 

7 7 ^ and -, = 2,3,... ,n- 1 . 58 

{Pj-lPj> iPjPn) [PjPn) 

We define the m-dimensional convex polytope 

Pm = {(Xl,... ,Xm) e [-1,1]" I \Xi-Xj\ < 1 Vi,j = 1, ... ,m} , (59) 

and replace the generation of the the ^-variables by the following: 

Algorithm 6 (IMPROVEMENT) 

1. generate {xi,X2, ■ ■ ■ ,x„^) distributed uniformly in P„^; 

2. define xo = and put, 
for alH = 2, . . . , n — 1. 



Because all the variables Xj are distributed uniformly such that \xi — Xj\ < 1, all quotients 
of ( |5^ ) are distributed such that they are between and ^m- In terms of the variables 
Xi, this means that we generate the volume of P^^, which is + 1, instead of the volume 
of [—1, 1]"«, which is 2"«. In |jl^, we give the algorithm to generate variables distributed 



uniformly in P^- We have to note here that this improvement only makes sense because 
the algorithm to generate these variables is very efficient. The total density changes such, 
that the function gn in Eq. ( ^OD has to be replaced by 

= (n^ + lKlog-U"^ "'"--""'''^''"'' • 

where the variables Xi are functions of the variables 2 as defined by (^). Because hit- 
and-miss MC is used to restrict generated events to cut PS, again only the normalization 
has to be calculated for the weight of an event. 

With this improvement, still a large number of events does not pass the cuts. The 
situation with PS is depicted in Fig.|I|. Phase space contains generated phase space which 
contains cut phase space. The problem is that most events fall in the shaded area, which 
is the piece of generated PS that is not contained in cut PS. To get a higher percentage of 
accepted events, we use a random variable G [0,^ni], instead of the fixed number ^m, to 
generate the variables ^[ 2- This means that the size of the generated PS becomes variable. 
If this is done with a probability distribution such that C,v can, in principle, become equal 
to C,^, then whole of cut phase space is still covered. We suggest the following, tunable, 
density: 

h (c \ an^ + l (log^^)""g ^ c ^ c \ ^ n (ao\ 

"^^'^ = (loge^)-e+i < < U , « > . (62) 
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phase space 

generated phase space 

cut phase space 



Figure 1: Schematic view on phase space. 

If a = 0, then log^t, is distributed uniformly in [0, log^m], and for larger a, the distribution 
peaks more and more towards C,v = ^m- Furthermore, the variable is easy to generate and 
the total generated density can be calculated exactly: gni^m] {(,}) should be replaced by 



^ ""^ + 1 (63) 



where ^low is the maximum of the ratios of scalar products in 



7 Results and conclusions 

We compare SARGE with RAMBO in the integration of the SPHEL-integrand for processes of 
the kind gg — >■ ng, which is given by 

y (64) 

(P1P2) (P2P3) (PsPi) ■ ■ ■ iPnPn+l) (Pn+lPn+2) (Pn+2Pl) 

where pi and p2 are the incoming momenta, and the first sum is over all permutations of 
(2, 3, . . . , n + 2) except the cyclic permutations. The results are presented in Tab.|^. The 
calculations were done at a CM-energy ^/s = 1000 with cuts pt = 40 on each transverse 
momentum and 6q = 30° on the angles between the momenta. We present the results for 
n = 4,5, 6, 7, calculated with RAMBO and SARGE with different values for a (Eq. (|6|)). The 
value of a is the estimate of the integral at an estimated error of 1% for n = 4,5,6 and 3% 



n 


4 


5 


6 


7 


T"sphel(s) 


5.40x10-5 


2.70x10-^ 


1.80x10-3 


1.41x10-2 


''exact (s) 


3.07x10"^ 


1.08 


3.35 


10.92 



Table 2: cpu-times (tsphel) in seconds needed to evaluate the SPHEL-integrand one time with a 300-MHz 
UhraSPARC-IIi processor, and the cpu-times (Tcxact) needed to evaluate the exact integrand, estimated 
with the help of Tab. 
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gg ^ 5g 
1% error 



gg 6g 
1% error 



gg 7g 
3% error 



alg. 


RAMBO 


SARGE, a = 0.0 


SARGE, a = 0.5 


SARGE, a = 10.0 


cr 


4.30 X 10^ 


4.31 X 10^ 


4.37x 10*^ 


4.32 X 10^ 




4, 736, 672 


296,050 


278, 702 


750,816 


iVac 


3,065,227 


111,320 


40,910 


23,373 




0.198 


0.0254 


0.0172 


0.0348 


icxa(h) 


262 


9.52 


3.51 


2.03 




RAMBO 


SARGE, a = 0.0 


SARGE, a; = 0.5 


SARGE, a = 10.0 




3.78 X 10^° 


3.81 X 10^° 


3.80 X 10^° 


3.81x10^° 




4, 243, 360 


715,585 


1,078, 129 


6,119,125 


iVac 


1,712,518 


167, 540 


36,385 


21,111 


^cpu (h) 


0.286 


0.133 


0.0758 


0.277 


^exa(h) 


514 


51.6 


11.7 


9.10 


als. 


RAMBO 


SARGE, a = 0.0 


SARGE, a = 0.5 


SARGE, a = 10.0 


cr 


3.07x 10^2 


3.05 X 10^2 


3.13x10^^ 


3.05 X 10^2 


gQ 


3,423,981 


2, 107, 743 


6, 136,375 


68,547,518 


iVac 


700, 482 


276, 344 


34,095 


17,973 


tcpu(h) 


0.685 


1.32 


0.471 


3.17 


icxa(h) 


653 


258 


32.2 


19.9 


alg. 


RAMBO 


SARGE, a = 0.0 


SARGE, a = 0.5 


SARGE, a = 10.0 


(J 


2.32 xlO^^ 


2.16x10^^ 


2.20 xlO^^ 


2.28 xlO^^ 




605,514 


710,602 


5,078,153 


125,471,887 


iVac 


49,915 


42, 394 


3,256 


1,789 


icpu(h) 


0.224 


1.86 


0.452 


6.74 


^exa(h) 


152 


130 


10.3 


12.2 



Table 3: Results for the integration of the SPHEL-integrand. The CM-energy and the cuts 
used are y/s = 1000, px = 40 and 6*0 = 30°. Presented are the finial result (cr), the 
number of generated (iVge) and accepted {Na,c) events, the cpu-time (tcpu) in hours, and 
the cpu-time (texa) it would take to integrate the exact matrix element, estimated with the 
help of Tab.|^. In the calculation of this table, adaptive multichanneling in the two cases 
of Section |5]^ was used, and S = 0.01 (Section |5.1|) . 
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4.4 
4.2 
4.0 
3.8 
3.6 
3.4 



T" 



T" 



T" 



T" 



1 1 T" 

RAMBO 

result integral (xlO^^° 




J I I I I I I L 



4.4 
4.2 
4.0 



1 



2 4 6 8 10 12 14 16 1^ 

iVac/lO^ 



3.4 



~~i 

SARGE 



result integral (xlO 



-io\ 



10 15 

iVac/lO^ 




6 8 10 12 14 16 18 

iVac/10^ 



10 15 

iVac/lO^ 



20 



25 




25 



Figure 2: The convergence of the MC-process in the integration of the SPHEL-integrand 
for n = 5, with y/s = 1000, pt = 40 and Oq = 30°. The upper graphs show the integral 
itself as function of the number of accepted events, together with the estimated bounds on 
the expected deviations. The lower graphs show the relative error. SARGE was used with 
adaptive multi-channeling in the two cases of Section |5.2| , with 6 = 0.01 (Section p.l| ) and 
without the variable C,v The number of generated events was 6, 699, 944, and the cpu-time 
was 0.308 hours. 
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for n = 7. These numbers are only printed to show that different results are compatible. 
Remember that they are not the whole cross sections: flux factors, color factors, sums and 
averages over helicities, and coupling constants are not included. The other data are the 
number of generated events (iVge), the number of accepted events (A^ac) that passed the 
cuts, the cpu-time consumed (tcpu), and the cpu-time the calculation would have consumed 
if the exact matrix element had been used (texa); both in hours. This final value is estimated 
with the help of Tab. § and the formula 

^exa = ^cpu + -^aclTcxact ~ TspHEl) 5 (65) 

where Xexact and Tsphel are the cpu-times it takes to evaluate the squared matrix element 
once. Remember that the integrand only has to be evaluated for accepted events. The 
calculations have been performed with a single 300-MHz UltraSPARC-IIi processor. 

The first conclusion we can draw is that SARGE outperforms RAMBO in computing time 
for all processes. This is especially striking for lower number of outgoing momenta, and 
this behavior has a simple explanation: we kept the CM-energy and the cuts fixed, so that 
there is less energy to distribute over the momenta if n is larger, and the cuts become 
relatively tighter. As a result, RAMBO gains on SARGE if n becomes larger. This effect would 
not appear if the energy, or the cuts, would scale with n like in Tab.[^. Another indication 
for this effect is the fact that the ratio Ng^^/Ngc for RAMBO, which estimates the ratio of the 
volumes of cut PS and whole PS, decreases with n. 

Another conclusion that can be drawn is that SARGE performs better if a is larger. 
Notice that the limit of a — oo is equivalent with dropping the improvement of the 
algorithm using the variable ^„ (Eq. (^)). Only if the integrand becomes too fiat, as in 
the case of n = 7 with the energy and the cuts as given in the table, smaller values are 
preferable. Then, too many events do not pass the cuts if a is large. 

As an extra illustration of the performance of SARGE, we present in Fig.H the evaluation 
of MC-integrals as function of the number of accepted events. Depicted are the integral a 
with the bounds on the expected deviation coming from the estimated expected error, and 
the relative error. Especially the graphs with the relative error are illustrative, since they 
show that it converges to zero more smoothly for SARGE then for RAMBO. Notice the spike 
for RAMBO around N^c = 25000, where an event obviously hits a singularity. 



8 Other pole structures 

The APS of (|I]) is not the only pole structure occurring in the squared amplitudes of QCD- 
processes; not even in purely gluonic processes. For example, in the case of gg 4g, also 
permutations of 

2 (66) 



{PlP3){P2P4)iP0Pl)iP0P2)iP0 - Pi -P2 
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occur 0. If one is able to generate momenta with this density, it can be included in the 
whole density with the use of the adaptive multichannel technique. In the interpretation 
of the transition amplitude as a sum of Feynman diagrams, this kind of pole structures 
typically come from t-channel diagrams, which are of the type 




and where, for this case, Qi = Pi + Ps and Q2 = P2 + Pa, so that k = pq — pi — p^. The 
natural way to generate a density with this pole structure is by generating Si = with a 
density proportional to l/sj, a variable t that plays the role of {po — pi — Ps)^, construct 
with this and some generated angles the momenta Qi, and then split new momenta from 
each of these. For n = 4, only two momenta have to split off each Qi, and there is a 
reasonable simple algorithm to generate these. 

We shall now just present the algorithm, and then show its correctness using the UAF. 
If we mention the generation of some random variable x 'with a density /(x)' in the 
following, we mean a density that is proportional to f{x), and we shall not always write 
down the normalization explicitly. Furthermore, s denotes the square of the CM-energy 
and A = A(s, Si, S2) the usual Mandelstam variable 

A = s"^ + sl + si — 2ssi — 2ss2 — 2siS2 ■ (67) 

Of course, a cut has to be implemented in order to generate momenta following (|66|), and 
we shall be able to put (piPj) > ^sq for the scalar products occurring in the denominator, 
where Sq only has to be larger than zero. To generate the momenta with density (|66D , one 
should 

Algorithm 7 (T-CHANNEL) 

1. generate si and S2 between sq and s with density 1/si and I/S2; 

2. generate t between s — Si — S2 ± a/ A(s, Si, S2) with density l/[t{t + 2si)(t + 2S2)]; 

3. put z ^ {s — Si — S2 — t)/VX and generate (p uniformly in [0, 27r); 

4. put Qi ^ (a/si + A/(4s), ^/\/JAs)n{z, if) ) and Q2 ^ ^/scq - Qi, 

5. for z = 1,2, generate Zi > 1 — 4so/(i + 2sj) with density 1/(1 — Zi) and cpi uniformly 
in [0,27r), and put qi ^ ly/s'i{l,n{zi,(fi) ); 

6. for i = 1,2, rotate to the CMF of Qi, then boost it to the CMF of Qi + Q2 to 
obtain pi, and put pi+2 Qi- Pi] 
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As a final step, the incoming momenta can be put to po ^ (cq + 63) and po <— 
|-\/s(eo — 63). The variables and Zi can easily be obtained by inversion (cf. 0). The 
variable t can best be obtained by generating x = log(2^siS2) — logt with the help of 
the rejection method (cf. [Q). In the UAF, the steps of the algorithm read as follows. 
Denoting 



£1 = Co + 63 



and 



nrm(s,Si,S2) = J 



62 = eo- 63 



dt 

t{t + 2si){t + 2s2] 



S — Si — S2 ± 



(6J 



we have 
1. 

2. 
3. 
4. 

5. 



1/4 

Sl - S2 



dSi ds2 0{so < Si^2 < s) 



llo 1 + 2^2//^- 1 l + 2si/h^ 
S2 1 + 252//;.+ Sl I + 2si/h^ 



(69) 



Sl S2 (log^) 
dt 



e{h_ <t<h+ 
t{t + 2si){t + 2S2) nrm(s, Si,S2) 
s — Sl — S2 — t\ dtp 

7x ; 2^ 



dz5 { z 



fi'QiMQi-\/^i 



Hi 

2 



dz, 0{l -Zi> 



4s 



4so 



5' [Q 

dipi 



h{z, Lp) ] d'^Q2 S'^iQi + Q2- \fs eo) 



log 



t+2s 



2so 



^ d% 5(g° - 1 Vi") 5^( - g°n(z., v^.) ) 



6. / n d'^b, 6\h - Hq^b,) 6\p, - nQ]n,^qi) 6\p,+2 + P^ - Qi) ■ 
^ i=i 

The various assignments imply the following identities. First of all, we have 

{Pi+Pi+2f = Q\ = Si . 

Using that 4ssi + A = (s + Si — ^2)^ we find 

(erQi) = s + Sl - S2 - z\f\ = t + 2si 
and the same for (1 *^2\ so that 

t = 4(po-Qi) - 2(pi + P3)' = -2{po -Pi- Psf . 
Denote £q. = TZb^HQ., so that = Cq-Pi. Because Cq-Ei ~ ei, we find that 

2{ei- qi 



1 - Zi 



(er^Q^Qi) (erQi) 



(70) 



(71) 



(72) 



(73) 



20 



so that 



{t + 2si){l-zi) = 8{po-pi) and {t + 2s2){l - Z2) = Sipo-P2) ■ (74) 

We can conclude so far that the algorithm generates the correct pole structure. For the 
further evaluation of the integrals one can forget about the factors Sj, t, t + 2sj and 1 — Zi 
in the denominators. Using that 

d% 5{q^ - iyi") 6-'{ I - q^h{z„ ^^)) ^ 2 d%^{qi) - n{z,, ^0 ) , (75) 

and replacing step 4 by 
2 



J] 2 Jsi + A d^Q.^sM) H4Qi) - ^) HviQi) - d\Qi + Q2 - v^eo) , (76) 



.1=1 



the integrals can easily be performed backwards, i.e., in the order gfj, Zj, foj, Qj, ^p, z, t, 
Si, S2- The density finally is 

Q ^^^y^^ g(2(poPi) > so)9{2{PoP2) > so)9{2ipm) > ^0)^(2(^2^4) > ^o) 



X 



iP0Pl){P0P2)iPlP3)iP2P4)[-iP0 -Pl~ P3, 

s f. s\\ t + 2si. t + 2s2 , . 
log — log — log — nrm(s, Si, S2) 

V So/ 2so 2so 



24(27r)3 

where = (p^ + Pi+2)^ and t = -2(po - Pi - Ps)^ 



(77) 



9 Appendices 

Appendix A 

We have to calculate the integral 

2.^f^V/da;ci^65(5^-l)^(60>0)^expf-£^60^ - 2r(2n)5(n) 



21:) j ^ ' ^ ' \ X ) (27r)"s'*-2 



where 

//•CK) 

The 'Euler substitution' &° = | (t'^^^ + v"^/^) casts the integral in the form 



/oo 
dv 



(^ + 1) 



2n 
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By the transformation v ^ 1/v it can easily be checked that the integral from 1 to oo is 
precisely equal to that from to 1, so that we may write 



2 7o + Vi2n 

where we have used, by writing z = 1/(1 + t>), that 



oo /•! 



JO 



References 



[1 

[2] 
[3] 

[4] 

[5] 

[6] 
[7] 



[9 
[10 

[11 



D.E. Knuth, The Art of Computer Programming, Vol.2. 2d ed. (Princeton, 1991). 

L. Devroye, Non-Uniform Random Variate Generation (Springer, 1986). 

R. Kleiss and R. Pittau, Weight optimization in multichannel Monte Carlo, Comp. 
Phys. Comm. 83 (1994) 141-146. 

W.J. Stirling, R. Kleiss and S.D. Ellis, A new Monte Carlo treatment of multiparticle 
phase space at high energy, Comp. Phys. Comm. 40 (1986) 359. 

J. CM. Kuijf, Multiparton production at hadron colliders, PhD thesis. University of 
Leiden, 1991. 

F. Berends and H. Kuijf, Jets at the LHC, Nucl. Phys. B353 (1991) 59-86. 

P. Draggiotis, R. Kleiss and C.G. Papadopoulos, On the computation of multigluon 
amplitudes, Nucl. Phys. B439 (1998) 157-164. 

F. Caravaglios, M.L. Mangano, M. Moretti and R. Pittau, A new approach to multi-jet 
calculations in hadron collisions, Nucl. Phys. B539 (1999) 215-232. 

A. van Hameren, R. Kleiss and P. Draggiotis, SARGE: an algorithm for generating 
QCD antennas, Phys. Lett. B 483 (2000) 124-130. 

A. van Hameren and R. Kleiss, A fast algorithm for generating a uniform distribution 
inside a high- dimensional polytope, preprint |physics/0003078 . 



P. Draggiotis, private communication. 



22 



